5.1.1 Prepare data

This section sets the variables to be used and prepares the data if necessary.

# load data 
load("R_objects/animal_data.Rdata")
params <- readRDS("R_objects/animal_params.RDS")
dat.clean <- dat

# Set names of variables
PREDICTOR <- "feedtreat"
OUTCOME <- "bw_gain08"
SUBJECT <- "rat_name"

# Create formula
PREDICTOR.F <- ifelse(length(PREDICTOR) > 1, paste(PREDICTOR, collapse = "*"), PREDICTOR)
FORMULA <- as.formula(paste(OUTCOME,PREDICTOR.F, sep = " ~ "))

# Summary samples in groups
dat.clean %>% group_by(across(all_of(PREDICTOR))) %>% get_summary_stats(!!sym(OUTCOME), type = c("mean_sd"))

## # A tibble: 4 × 5
##   feedtreat variable      n  mean    sd
##   <chr>     <fct>     <dbl> <dbl> <dbl>
## 1 HF_CTRL   bw_gain08    12  21.3  3.46
## 2 HF_PFOS   bw_gain08    12  19.4  3.46
## 3 LF_CTRL   bw_gain08    12  24.1  3.54
## 4 LF_PFOS   bw_gain08    12  23.5  3.08

dat.clean %>% group_by(across(all_of(PREDICTOR))) %>% get_summary_stats(!!sym(OUTCOME), type = c("full"))

## # A tibble: 4 × 14
##   feedtreat variable      n   min   max median    q1    q3   iqr   mad  mean
##   <chr>     <fct>     <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 HF_CTRL   bw_gain08    12  16.9  26.7   20.2  18.6  24.3  5.69  3.43  21.3
## 2 HF_PFOS   bw_gain08    12  12.8  24.5   19.6  17.8  21.5  3.7   2.96  19.4
## 3 LF_CTRL   bw_gain08    12  18.4  28.7   23.9  21.0  27.2  6.17  4.65  24.1
## 4 LF_PFOS   bw_gain08    12  18.7  27.7   23.7  21.7  25.8  4.09  3.55  23.5
## # ℹ 3 more variables: sd <dbl>, se <dbl>, ci <dbl>

# Create plot
bxp <- dat.clean %>%
  ggboxplot(x = if_else(length(PREDICTOR) > 1, PREDICTOR[2],PREDICTOR[1]),
            y = OUTCOME,
            color = PREDICTOR[1],
            facet.by = if(length(PREDICTOR) == 3) PREDICTOR[3],
            palette = params$COL)
bxp

Assumptions and preliminary tests

The ANOVA tests assume the following characteristics about the data:

• Independence of the observations. Each subject should belong to only one group. There is no relationship between the observations in each group.
  This is already done for the whole project

• No significant outliers in the two groups

• Normality. the data for each group should be approximately normally distributed.

• Homogeneity of variances. the variance of the outcome variable should be equal in each group.

In this section, we’ll perform some preliminary tests to check whether these assumptions are met.

Identify outliers
Outliers can be easily identified using boxplot methods, implemented in the R function identify_outliers() [rstatix package].

# Test for outliers
dat.clean %>% 
  group_by(across(all_of(PREDICTOR))) %>% 
  identify_outliers(!!sym(OUTCOME))

##  [1] feedtreat         rat_org           rat_name          dissection       
##  [5] cage              treatment         feed              feedtreatday     
##  [9] bw_minus18        bw_0              bw_1              bw_2             
## [13] bw_3              bw_4              bw_5              bw_6             
## [17] bw_7              bw_8              bw_12             bw_16            
## [21] bw_21             bw_gain08         bw_gain021        bw_gain821       
## [25] bloodvol_8        bloodvol_16       bloodvol_21       cecum_wt         
## [29] cecum_wtbw        cecum_norm        liver_wt          liver_wtbw       
## [33] liver_norm        brain_wt_estimate transit_0         transit_7        
## [37] transit_20        dose_total_ml     pfos_total_mg     pfos_liver_ugg   
## [41] pfos_liver_mg     pfos_liver8_pct   pfos_brain_ugg    pfos_brain_uggbw 
## [45] pfos_brain_ug     pfos_brain_pct    pfos_serum8_ugml  pfos_serum8_mg   
## [49] pfos_serum8_pct   pfos_serum16_ugml pfos_serum16_mg   pfos_serum16_pct 
## [53] pfos_serum21_ugml pfos_serum21_mg   pfos_serum21_pct  pfos_feces8_ugg  
## [57] pfos_feces12_ugg  pfos_feces16_ugg  pfos_feces21_ugg  pfos_cecum_ugg   
## [61] pfos_urine8_ugml  pfos_urine16_ugml pfos_urine21_ugml pH_je_up         
## [65] pH_je_down        pH_ileum          pH_cecum          acetic           
## [69] formic            propanoic         m2_propanoic      butanoic         
## [73] m3_butanoic       pentanoic         m4_pentanoic      hexanoic         
## [77] heptanoic         is.outlier        is.extreme       
## <0 rækker> (eller 0-længde row.names)

Data contains two outliers: sample from rat_name R01 and R30.

Check normality
QQ plot and Shapiro-Wilk test of normality are used to analyze the model residuals.

# Build the linear model
model  <- lm(FORMULA, data = dat.clean)
# Create a QQ plot of residuals
ggqqplot(residuals(model))

# Compute Shapiro-Wilk test of normality
shapiro_test(residuals(model))

## # A tibble: 1 × 3
##   variable         statistic p.value
##   <chr>                <dbl>   <dbl>
## 1 residuals(model)     0.961   0.114

Test homogneity of variance assumption
1. The residuals versus fits plot can be used to check the homogeneity of variances.

plot(model, 1)

2. It’s also possible to use the Levene’s test to check the homogeneity of variances:

dat.clean %>% levene_test(FORMULA)

## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.

## # A tibble: 1 × 4
##     df1   df2 statistic     p
##   <int> <int>     <dbl> <dbl>
## 1     3    44     0.133 0.940

# Save result
EQUAL.VAR <- dat.clean %>% levene_test(FORMULA) %>% pull(p) > 0.05

## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.

This shows that body weight gain data has two outliers, has equal variance and is normally distributed without the outliers according to Shapiro-Wilk test. Therefore we use a one-way ANOVA test with Tukey’s honest significance test.

5.1.2 ANOVA One-Way test

Perform test

If we had equality of variance we can now run a one-way ANOVA tests anova_test() (if we have equal variance) or a welch_anova_test() (if variance vary).

if(EQUAL.VAR) {
  res.aov <- dat.clean %>% anova_test(FORMULA)
  res.aov
} else {
  res.aov <- dat.clean %>% welch_anova_test(FORMULA)
  res.aov
}

## ANOVA Table (type II tests)
## 
##      Effect DFn DFd     F     p p<.05   ges
## 1 feedtreat   3  44 4.863 0.005     * 0.249

Perform posthoc test

A significant one-way ANOVA is generally followed up by Tukey post-hoc tests to perform multiple pairwise comparisons between groups. When running relaxed Welch one-way test, the Games-Howell post hoc test or pairwise t-tests (with no assumption of equal variances) can be used to compare all possible combinations of group differences.

if(EQUAL.VAR) {
  pwc <- dat.clean %>% tukey_hsd(FORMULA)
  pwc
} else {
  pwc <- dat.clean %>% games_howell_test(FORMULA)
  pwc
}

## # A tibble: 6 × 9
##   term      group1  group2  null.value estimate conf.low conf.high   p.adj
## * <chr>     <chr>   <chr>        <dbl>    <dbl>    <dbl>     <dbl>   <dbl>
## 1 feedtreat HF_CTRL HF_PFOS          0   -1.88    -5.57       1.82 0.533  
## 2 feedtreat HF_CTRL LF_CTRL          0    2.84    -0.851      6.54 0.184  
## 3 feedtreat HF_CTRL LF_PFOS          0    2.21    -1.49       5.91 0.391  
## 4 feedtreat HF_PFOS LF_CTRL          0    4.72     1.03       8.42 0.00737
## 5 feedtreat HF_PFOS LF_PFOS          0    4.09     0.392      7.78 0.025  
## 6 feedtreat LF_CTRL LF_PFOS          0   -0.634   -4.33       3.06 0.968  
## # ℹ 1 more variable: p.adj.signif <chr>

5.1.3 Create nested figure

# Set variables for inner and outer analysis, and variable for facet
INNER.VAR <- "treatment"
OUTER.VAR <- "feed"
FACETVAR <- "dissection"

# Statistics costumed for facet plotting
## Pairwise comparison for inner variable
stat.in <- dat.clean %>%
  group_by(.data[[OUTER.VAR]]) %>%
  t_test(as.formula(paste(OUTCOME,"~",INNER.VAR, sep = " ")), 
         paired = FALSE, var.equal = EQUAL.VAR, 
         detailed = TRUE, alternative = "two.sided") %>%
  # adjust_pvalue(method = "bonferroni") %>%
  add_significance() %>%
  p_format("p", accuracy = 0.001, trailing.zero = TRUE, new.col = TRUE)
stat.in

## # A tibble: 2 × 18
##   feed  estimate estimate1 estimate2 .y.     group1 group2    n1    n2 statistic
## * <chr>    <dbl>     <dbl>     <dbl> <chr>   <chr>  <chr>  <int> <int>     <dbl>
## 1 HF       1.88       21.3      19.4 bw_gai… CTRL   PFOS      12    12     1.33 
## 2 LF       0.634      24.1      23.5 bw_gai… CTRL   PFOS      12    12     0.468
## # ℹ 8 more variables: p <dbl>, df <dbl>, conf.low <dbl>, conf.high <dbl>,
## #   method <chr>, alternative <chr>, p.signif <chr>, p.format <chr>

## Pairwise comparison for outer variable
stat.out <- dat.clean %>%
  # group_by(.data[[FACETVAR]]) %>%
  t_test(as.formula(paste(OUTCOME,"~", OUTER.VAR, sep = " ")), 
         paired = FALSE, var.equal = EQUAL.VAR, 
         detailed = TRUE, alternative = "two.sided") %>%
  # adjust_pvalue(method = "bonferroni") %>%
  add_significance() %>%
  p_format("p", accuracy = 0.001, trailing.zero = TRUE, new.col = TRUE)
stat.out

## # A tibble: 1 × 17
##   estimate estimate1 estimate2 .y.    group1 group2    n1    n2 statistic      p
## *    <dbl>     <dbl>     <dbl> <chr>  <chr>  <chr>  <int> <int>     <dbl>  <dbl>
## 1    -3.47      20.3      23.8 bw_ga… HF     LF        24    24     -3.54 9.3e-4
## # ℹ 7 more variables: df <dbl>, conf.low <dbl>, conf.high <dbl>, method <chr>,
## #   alternative <chr>, p.signif <chr>, p.format <chr>

## Calculate positions statistics on plot
stat.in <- stat.in %>% add_xy_position(x = OUTER.VAR, dodge = 0.8)
stat.out <- stat.out %>% add_xy_position(x = OUTER.VAR)
stat.out$y.position <- max(stat.in$y.position)*1.1

# Sort dat.clean
dat.clean <- dat.clean[order(dat.clean$feedtreat),]

# Create plot
p <- ggboxplot(dat.clean, x = OUTER.VAR, y = OUTCOME,
               color = PREDICTOR,
               fill = PREDICTOR, 
               add =  "jitter", 
               add.params = list(size = 1)) +
  scale_fill_manual(values = params$COL1, name = "Group", labels = c("HF_CTRL" = "HF-CTRL","HF_PFOS" = "HF-PFOS","LF_CTRL" = "LF-CTRL","LF_PFOS" = "LF-PFOS")) +
  scale_color_manual(values = c("black","black","black","black")) +
  scale_y_continuous(name = "Bodyweight gain",limits = c(10,35),breaks = seq(10,35,5), labels = function(x) paste0(x, "%")) +
  guides(color = "none") +
  theme(axis.title.x = element_blank())

p <- p + stat_pvalue_manual(stat.in, label = "p.format",tip.length = 0, hide.ns = FALSE)+
  stat_pvalue_manual(stat.out, label = "p.format", tip.length = 0, hide.ns = FALSE)
p

# Plot for saving without legend
p2 <- p + theme(legend.position = "none")

# Output plot
ggsave(filename = paste0("plots/animal_data/weight/",OUTCOME,"_plot.png"), p2, device = "png", dpi = 300, units = "mm", width = 100, height = 100)
ggsave(filename = paste0("plots/animal_data/weight/",OUTCOME,"_plot.pdf"), p2, device = "pdf", dpi = 300, units = "mm", width = 100, height = 100)

# clear the environment and release memory
rm(list = ls(all.names = TRUE)) #will clear all objects includes hidden objects.
invisible(gc()) #free up memory and report the memory usage.

5.2.1 Prepare data

This section sets the variables to be used and prepares the data if necessary.

# load data 
load("R_objects/animal_data.Rdata")
params <- readRDS("R_objects/animal_params.RDS")

# Subset
dat.clean <- subset(dat, !is.na(bw_gain021))

# Set names of variables
PREDICTOR <- "feedtreat"
OUTCOME <- "bw_gain021"
SUBJECT <- "rat_name"

# Create formula
PREDICTOR.F <- ifelse(length(PREDICTOR) > 1, paste(PREDICTOR, collapse = "*"), PREDICTOR)
FORMULA <- as.formula(paste(OUTCOME,PREDICTOR.F, sep = " ~ "))

# Summary samples in groups
dat.clean %>% group_by(across(all_of(PREDICTOR))) %>% get_summary_stats(!!sym(OUTCOME), type = "mean_sd")

## # A tibble: 4 × 5
##   feedtreat variable       n  mean    sd
##   <chr>     <fct>      <dbl> <dbl> <dbl>
## 1 HF_CTRL   bw_gain021     6  45.5  6.86
## 2 HF_PFOS   bw_gain021     6  47.4  5.93
## 3 LF_CTRL   bw_gain021     6  57.0  7.82
## 4 LF_PFOS   bw_gain021     6  53.0  5.48

dat.clean %>% group_by(across(all_of(PREDICTOR))) %>% get_summary_stats(!!sym(OUTCOME), type = c("full"))

## # A tibble: 4 × 14
##   feedtreat variable       n   min   max median    q1    q3   iqr   mad  mean
##   <chr>     <fct>      <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 HF_CTRL   bw_gain021     6  36.5  54.6   46.8  40.3  49.3  8.92  7.71  45.5
## 2 HF_PFOS   bw_gain021     6  36.6  52.1   49.5  45.7  51.3  5.58  3.67  47.4
## 3 LF_CTRL   bw_gain021     6  48.9  68.1   54.3  51.6  62.7 11.1   6.40  57.0
## 4 LF_PFOS   bw_gain021     6  47.0  60.8   53.1  48.2  56.2  8.03  7.24  53.0
## # ℹ 3 more variables: sd <dbl>, se <dbl>, ci <dbl>

# Create plot
bxp <- dat.clean %>%
  ggboxplot(x = if_else(length(PREDICTOR) > 1, PREDICTOR[2],PREDICTOR[1]),
            y = OUTCOME,
            color = PREDICTOR[1],
            facet.by = if(length(PREDICTOR) == 3) PREDICTOR[3],
            palette = params$COL)
bxp

# Test for outliers
dat.clean %>% 
  group_by(across(all_of(PREDICTOR))) %>% 
  identify_outliers(!!sym(OUTCOME))

## # A tibble: 1 × 79
##   feedtreat rat_org rat_name dissection  cage treatment feed  feedtreatday
##   <chr>     <chr>   <chr>    <chr>      <int> <chr>     <chr> <chr>       
## 1 HF_PFOS   R46     R46      d21           23 PFOS      HF    HF_PFOS_d21 
## # ℹ 71 more variables: bw_minus18 <dbl>, bw_0 <dbl>, bw_1 <int>, bw_2 <int>,
## #   bw_3 <int>, bw_4 <int>, bw_5 <int>, bw_6 <int>, bw_7 <int>, bw_8 <int>,
## #   bw_12 <dbl>, bw_16 <dbl>, bw_21 <dbl>, bw_gain08 <dbl>, bw_gain021 <dbl>,
## #   bw_gain821 <dbl>, bloodvol_8 <dbl>, bloodvol_16 <dbl>, bloodvol_21 <dbl>,
## #   cecum_wt <dbl>, cecum_wtbw <dbl>, cecum_norm <dbl>, liver_wt <dbl>,
## #   liver_wtbw <dbl>, liver_norm <dbl>, brain_wt_estimate <dbl>,
## #   transit_0 <int>, transit_7 <int>, transit_20 <int>, dose_total_ml <dbl>, …

# Build the linear model
model  <- lm(FORMULA, data = dat.clean)
# Create a QQ plot of residuals
ggqqplot(residuals(model))

# Compute Shapiro-Wilk test of normality
shapiro_test(residuals(model))

## # A tibble: 1 × 3
##   variable         statistic p.value
##   <chr>                <dbl>   <dbl>
## 1 residuals(model)     0.969   0.636

plot(model, 1)

dat.clean %>% levene_test(FORMULA)

## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.

## # A tibble: 1 × 4
##     df1   df2 statistic     p
##   <int> <int>     <dbl> <dbl>
## 1     3    20     0.302 0.824

# Save result
EQUAL.VAR <- dat.clean %>% levene_test(FORMULA) %>% pull(p) > 0.05

## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.

This shows that body weight gain data has two outliers, has equal variance and is normally distributed without the outliers according to Shapiro-Wilk test. Therefore we use a one-way ANOVA test with Tukey’s honest significance test.

5.2.2 ANOVA One-Way test

Perform test

If we had equality of variance we can now run a one-way ANOVA tests anova_test() (if we have equal variance) or a welch_anova_test() (if variance vary).

if(EQUAL.VAR) {
  res.aov <- dat.clean %>% anova_test(FORMULA)
  res.aov
} else {
  res.aov <- dat.clean %>% welch_anova_test(FORMULA)
  res.aov
}

## ANOVA Table (type II tests)
## 
##      Effect DFn DFd     F     p p<.05   ges
## 1 feedtreat   3  20 3.817 0.026     * 0.364

Perform posthoc test

A significant one-way ANOVA is generally followed up by Tukey post-hoc tests to perform multiple pairwise comparisons between groups. When running relaxed Welch one-way test, the Games-Howell post hoc test or pairwise t-tests (with no assumption of equal variances) can be used to compare all possible combinations of group differences.

if(EQUAL.VAR) {
  pwc <- dat.clean %>% tukey_hsd(FORMULA)
  pwc
} else {
  pwc <- dat.clean %>% games_howell_test(FORMULA)
  pwc
}

## # A tibble: 6 × 9
##   term  group1 group2 null.value estimate conf.low conf.high  p.adj p.adj.signif
## * <chr> <chr>  <chr>       <dbl>    <dbl>    <dbl>     <dbl>  <dbl> <chr>       
## 1 feed… HF_CT… HF_PF…          0     1.86   -8.78      12.5  0.961  ns          
## 2 feed… HF_CT… LF_CT…          0    11.5     0.838     22.1  0.0316 *           
## 3 feed… HF_CT… LF_PF…          0     7.46   -3.18      18.1  0.235  ns          
## 4 feed… HF_PF… LF_CT…          0     9.62   -1.02      20.3  0.0853 ns          
## 5 feed… HF_PF… LF_PF…          0     5.61   -5.04      16.2  0.471  ns          
## 6 feed… LF_CT… LF_PF…          0    -4.02  -14.7        6.63 0.719  ns

5.2.3 Create nested figure

# Set variables for inner and outer analysis, and variable for facet
INNER.VAR <- "treatment"
OUTER.VAR <- "feed"
FACETVAR <- "dissection"

# Statistics costumed for facet plotting
## Pairwise comparison for inner variable
stat.in <- dat.clean %>%
  group_by(.data[[OUTER.VAR]]) %>%
  t_test(as.formula(paste(OUTCOME,"~",INNER.VAR, sep = " ")), 
         paired = FALSE, var.equal = EQUAL.VAR, 
         detailed = TRUE, alternative = "two.sided") %>%
  # adjust_pvalue(method = "bonferroni") %>%
  add_significance() %>%
  p_format("p", accuracy = 0.001, trailing.zero = TRUE, new.col = TRUE)
stat.in

## # A tibble: 2 × 18
##   feed  estimate estimate1 estimate2 .y.     group1 group2    n1    n2 statistic
## * <chr>    <dbl>     <dbl>     <dbl> <chr>   <chr>  <chr>  <int> <int>     <dbl>
## 1 HF       -1.86      45.5      47.4 bw_gai… CTRL   PFOS       6     6    -0.502
## 2 LF        4.02      57.0      53.0 bw_gai… CTRL   PFOS       6     6     1.03 
## # ℹ 8 more variables: p <dbl>, df <dbl>, conf.low <dbl>, conf.high <dbl>,
## #   method <chr>, alternative <chr>, p.signif <chr>, p.format <chr>

## Pairwise comparison for outer variable
stat.out <- dat.clean %>%
  # group_by(.data[[FACETVAR]]) %>%
  t_test(as.formula(paste(OUTCOME,"~", OUTER.VAR, sep = " ")), 
         paired = FALSE, var.equal = EQUAL.VAR, 
         detailed = TRUE, alternative = "two.sided") %>%
  # adjust_pvalue(method = "bonferroni") %>%
  add_significance() %>%
  p_format("p", accuracy = 0.001, trailing.zero = TRUE, new.col = TRUE)
stat.out

## # A tibble: 1 × 17
##   estimate estimate1 estimate2 .y.   group1 group2    n1    n2 statistic       p
## *    <dbl>     <dbl>     <dbl> <chr> <chr>  <chr>  <int> <int>     <dbl>   <dbl>
## 1    -8.54      46.4      55.0 bw_g… HF     LF        12    12     -3.23 0.00389
## # ℹ 7 more variables: df <dbl>, conf.low <dbl>, conf.high <dbl>, method <chr>,
## #   alternative <chr>, p.signif <chr>, p.format <chr>

## Calculate positions statistics on plot
stat.in <- stat.in %>% add_xy_position(x = OUTER.VAR, dodge = 0.8)
stat.out <- stat.out %>% add_xy_position(x = OUTER.VAR)
stat.out$y.position <- max(stat.in$y.position)*1.1

# Sort dat.clean
dat.clean <- dat.clean[order(dat.clean$feedtreat),]

# Create plot
p <- ggboxplot(dat.clean, x = OUTER.VAR, y = OUTCOME,
               color = PREDICTOR,
               fill = PREDICTOR, 
               add =  "jitter", 
               add.params = list(size = 1)) +
  scale_fill_manual(values = params$COL1, name = "Group", labels = c("HF_CTRL" = "HF-CTRL","HF_PFOS" = "HF-PFOS","LF_CTRL" = "LF-CTRL","LF_PFOS" = "LF-PFOS")) +
  scale_color_manual(values = c("black","black","black","black")) +
  scale_y_continuous(name = "Bodyweight gain",limits = c(30,80),breaks = seq(30,80,10), labels = function(x) paste0(x, "%")) +
  guides(color = "none") +
  theme(axis.title.x = element_blank())

p <- p + stat_pvalue_manual(stat.in, label = "p.format",tip.length = 0, hide.ns = FALSE)+
  stat_pvalue_manual(stat.out, label = "p.format", tip.length = 0, hide.ns = FALSE)

p

# Plot for saving without legend
p2 <- p + theme(legend.position = "none")

# Output plot
ggsave(filename = paste0("plots/animal_data/weight/",OUTCOME,"_plot.png"), p2, device = "png", dpi = 300, units = "mm", width = 100, height = 100)
ggsave(filename = paste0("plots/animal_data/weight/",OUTCOME,"_plot.pdf"), p2, device = "pdf", dpi = 300, units = "mm", width = 100, height = 100)

# clear the environment and release memory
rm(list = ls(all.names = TRUE)) #will clear all objects includes hidden objects.
invisible(gc()) #free up memory and report the memory usage.

5.3.1 Prepare data

This section sets the variables to be used and prepares the data if necessary.

# load data 
load("R_objects/animal_data.Rdata")
params <- readRDS("R_objects/animal_params.RDS")

# Subset
dat.clean <- subset(dat, !is.na(bw_gain821))

# Set names of variables
PREDICTOR <- "feedtreat"
OUTCOME <- "bw_gain821"
SUBJECT <- "rat_name"

# Create formula
PREDICTOR.F <- ifelse(length(PREDICTOR) > 1, paste(PREDICTOR, collapse = "*"), PREDICTOR)
FORMULA <- as.formula(paste(OUTCOME,PREDICTOR.F, sep = " ~ "))

# Summary samples in groups
dat.clean %>% group_by(across(all_of(PREDICTOR))) %>% get_summary_stats(!!sym(OUTCOME), type = "mean_sd")

## # A tibble: 4 × 5
##   feedtreat variable       n  mean    sd
##   <chr>     <fct>      <dbl> <dbl> <dbl>
## 1 HF_CTRL   bw_gain821     6  19.5  4.42
## 2 HF_PFOS   bw_gain821     6  22.0  5.08
## 3 LF_CTRL   bw_gain821     6  25.2  4.55
## 4 LF_PFOS   bw_gain821     6  23.9  4.47

# Create plot
bxp <- dat.clean %>%
  ggboxplot(x = if_else(length(PREDICTOR) > 1, PREDICTOR[2],PREDICTOR[1]),
            y = OUTCOME,
            color = PREDICTOR[1],
            facet.by = if(length(PREDICTOR) == 3) PREDICTOR[3],
            palette = params$COL)
bxp

# Test for outliers
dat.clean %>% 
  group_by(across(all_of(PREDICTOR))) %>% 
  identify_outliers(!!sym(OUTCOME))

## # A tibble: 4 × 79
##   feedtreat rat_org rat_name dissection  cage treatment feed  feedtreatday
##   <chr>     <chr>   <chr>    <chr>      <int> <chr>     <chr> <chr>       
## 1 HF_CTRL   R36     R36      d21           18 CTRL      HF    HF_CTRL_d21 
## 2 HF_PFOS   R46     R46      d21           23 PFOS      HF    HF_PFOS_d21 
## 3 LF_PFOS   R19     R19      d21           10 PFOS      LF    LF_PFOS_d21 
## 4 LF_PFOS   R23     R23      d21           12 PFOS      LF    LF_PFOS_d21 
## # ℹ 71 more variables: bw_minus18 <dbl>, bw_0 <dbl>, bw_1 <int>, bw_2 <int>,
## #   bw_3 <int>, bw_4 <int>, bw_5 <int>, bw_6 <int>, bw_7 <int>, bw_8 <int>,
## #   bw_12 <dbl>, bw_16 <dbl>, bw_21 <dbl>, bw_gain08 <dbl>, bw_gain021 <dbl>,
## #   bw_gain821 <dbl>, bloodvol_8 <dbl>, bloodvol_16 <dbl>, bloodvol_21 <dbl>,
## #   cecum_wt <dbl>, cecum_wtbw <dbl>, cecum_norm <dbl>, liver_wt <dbl>,
## #   liver_wtbw <dbl>, liver_norm <dbl>, brain_wt_estimate <dbl>,
## #   transit_0 <int>, transit_7 <int>, transit_20 <int>, dose_total_ml <dbl>, …

# Build the linear model
model  <- lm(FORMULA, data = dat.clean)
# Create a QQ plot of residuals
ggqqplot(residuals(model))

# Compute Shapiro-Wilk test of normality
shapiro_test(residuals(model))

## # A tibble: 1 × 3
##   variable         statistic p.value
##   <chr>                <dbl>   <dbl>
## 1 residuals(model)     0.955   0.348

plot(model, 1)

dat.clean %>% levene_test(FORMULA)

## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.

## # A tibble: 1 × 4
##     df1   df2 statistic     p
##   <int> <int>     <dbl> <dbl>
## 1     3    20    0.0776 0.971

# Save result
EQUAL.VAR <- dat.clean %>% levene_test(FORMULA) %>% pull(p) > 0.05

## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.

This shows that body weight gain data has two outliers, has equal variance and is normally distributed without the outliers according to Shapiro-Wilk test. Therefore we use a one-way ANOVA test with Tukey’s honest significance test.

5.3.2 ANOVA One-Way test

Perform test

If we had equality of variance we can now run a one-way ANOVA tests anova_test() (if we have equal variance) or a welch_anova_test() (if variance vary).

if(EQUAL.VAR) {
  res.aov <- dat.clean %>% anova_test(FORMULA)
  res.aov
} else {
  res.aov <- dat.clean %>% welch_anova_test(FORMULA)
  res.aov
}

## ANOVA Table (type II tests)
## 
##      Effect DFn DFd     F     p p<.05   ges
## 1 feedtreat   3  20 1.716 0.196       0.205

Perform posthoc test

A significant one-way ANOVA is generally followed up by Tukey post-hoc tests to perform multiple pairwise comparisons between groups. When running relaxed Welch one-way test, the Games-Howell post hoc test or pairwise t-tests (with no assumption of equal variances) can be used to compare all possible combinations of group differences.

if(EQUAL.VAR) {
  pwc <- dat.clean %>% tukey_hsd(FORMULA)
  pwc
} else {
  pwc <- dat.clean %>% games_howell_test(FORMULA)
  pwc
}

## # A tibble: 6 × 9
##   term   group1 group2 null.value estimate conf.low conf.high p.adj p.adj.signif
## * <chr>  <chr>  <chr>       <dbl>    <dbl>    <dbl>     <dbl> <dbl> <chr>       
## 1 feedt… HF_CT… HF_PF…          0     2.51    -4.99     10.0  0.786 ns          
## 2 feedt… HF_CT… LF_CT…          0     5.71    -1.79     13.2  0.177 ns          
## 3 feedt… HF_CT… LF_PF…          0     4.40    -3.09     11.9  0.378 ns          
## 4 feedt… HF_PF… LF_CT…          0     3.20    -4.29     10.7  0.637 ns          
## 5 feedt… HF_PF… LF_PF…          0     1.89    -5.60      9.39 0.893 ns          
## 6 feedt… LF_CT… LF_PF…          0    -1.31    -8.80      6.19 0.961 ns

5.3.3 Create figure

# Set variables for inner and outer analysis, and variable for facet
INNER.VAR <- "treatment"
OUTER.VAR <- "feed"
FACETVAR <- "dissection"

# Statistics costumed for facet plotting
## Pairwise comparison for inner variable
stat.in <- dat.clean %>%
  group_by(.data[[OUTER.VAR]]) %>%
  t_test(as.formula(paste(OUTCOME,"~",INNER.VAR, sep = " ")), 
         paired = FALSE, var.equal = EQUAL.VAR, 
         detailed = TRUE, alternative = "two.sided") %>%
  # adjust_pvalue(method = "bonferroni") %>%
  add_significance() %>%
  p_format("p", accuracy = 0.001, trailing.zero = TRUE, new.col = TRUE)
stat.in

## # A tibble: 2 × 18
##   feed  estimate estimate1 estimate2 .y.     group1 group2    n1    n2 statistic
## * <chr>    <dbl>     <dbl>     <dbl> <chr>   <chr>  <chr>  <int> <int>     <dbl>
## 1 HF       -2.51      19.5      22.0 bw_gai… CTRL   PFOS       6     6    -0.912
## 2 LF        1.31      25.2      23.9 bw_gai… CTRL   PFOS       6     6     0.502
## # ℹ 8 more variables: p <dbl>, df <dbl>, conf.low <dbl>, conf.high <dbl>,
## #   method <chr>, alternative <chr>, p.signif <chr>, p.format <chr>

## Pairwise comparison for outer variable
stat.out <- dat.clean %>%
  # group_by(.data[[FACETVAR]]) %>%
  t_test(as.formula(paste(OUTCOME,"~", OUTER.VAR, sep = " ")), 
         paired = FALSE, var.equal = EQUAL.VAR, 
         detailed = TRUE, alternative = "two.sided") %>%
  # adjust_pvalue(method = "bonferroni") %>%
  add_significance() %>%
  p_format("p", accuracy = 0.001, trailing.zero = TRUE, new.col = TRUE)
stat.out

## # A tibble: 1 × 17
##   estimate estimate1 estimate2 .y.    group1 group2    n1    n2 statistic      p
## *    <dbl>     <dbl>     <dbl> <chr>  <chr>  <chr>  <int> <int>     <dbl>  <dbl>
## 1    -3.80      20.8      24.6 bw_ga… HF     LF        12    12     -2.05 0.0525
## # ℹ 7 more variables: df <dbl>, conf.low <dbl>, conf.high <dbl>, method <chr>,
## #   alternative <chr>, p.signif <chr>, p.format <chr>

## Calculate positions statistics on plot
stat.in <- stat.in %>% add_xy_position(x = OUTER.VAR, dodge = 0.8)
stat.out <- stat.out %>% add_xy_position(x = OUTER.VAR)
stat.out$y.position <- max(stat.in$y.position)*1.1

# Sort dat.clean
dat.clean <- dat.clean[order(dat.clean$feedtreat),]

# Create plot
p <- ggboxplot(dat.clean, x = OUTER.VAR, y = OUTCOME,
               color = PREDICTOR,
               fill = PREDICTOR, 
               add =  "jitter", 
               add.params = list(size = 1)) +
  scale_fill_manual(values = params$COL1, name = "Group", labels = c("HF_CTRL" = "HF-CTRL","HF_PFOS" = "HF-PFOS","LF_CTRL" = "LF-CTRL","LF_PFOS" = "LF-PFOS")) +
  scale_color_manual(values = c("black","black","black","black")) +
  scale_y_continuous(name = "Bodyweight gain",limits = c(10,40),breaks = seq(10,40,5), labels = function(x) paste0(x, "%")) +
  guides(color = "none") +
  theme(axis.title.x = element_blank())

p <- p + stat_pvalue_manual(stat.in, label = "p.format",tip.length = 0, hide.ns = FALSE)+
  stat_pvalue_manual(stat.out, label = "p.format", tip.length = 0, hide.ns = FALSE)

p

# Plot for saving without legend
p2 <- p + theme(legend.position = "none")

# Output plot
ggsave(filename = paste0("plots/animal_data/weight/",OUTCOME,"_plot.png"), p2, device = "png", dpi = 300, units = "mm", width = 100, height = 100)
ggsave(filename = paste0("plots/animal_data/weight/",OUTCOME,"_plot.pdf"), p2, device = "pdf", dpi = 300, units = "mm", width = 100, height = 100)

# clear the environment and release memory
rm(list = ls(all.names = TRUE)) #will clear all objects includes hidden objects.
invisible(gc()) #free up memory and report the memory usage.

5.4.1 Prepare data and growth rates

# Load data
load("R_objects/animal_data.Rdata")
params <- readRDS("R_objects/animal_params.RDS")

# Set parameters
PREDICTOR <- c("bw_time_n","treatment","feed")
OUTCOME <- "wt_n"
SUBJECT <- "rat_name"

# Create data frame for data representation
dat.clean <- dat %>% select(rat_name, feed, treatment, feedtreat, dissection, bw_0, bw_1, bw_2, bw_3, bw_4, bw_5, bw_6, bw_7, bw_8, bw_12, bw_16, bw_21)

# Normalise data
dat.clean <- dat.clean %>% mutate(bw_0_n = bw_0 / bw_0, bw_1_n = bw_1 / bw_0, bw_2_n = bw_2 / bw_0, bw_3_n = bw_3 / bw_0, bw_4_n = bw_4 / bw_0, bw_5_n = bw_5 / bw_0, bw_6_n = bw_6 / bw_0, bw_7_n = bw_7 / bw_0, bw_8_n = bw_8 / bw_0, bw_12_n = bw_12 / bw_0, bw_16_n = bw_16 / bw_0, bw_21_n = bw_21 / bw_0)

# Pivot longer
dat.bwn <- dat.clean %>%
    pivot_longer(., cols = c(bw_0_n, bw_1_n, bw_2_n, bw_3_n, bw_4_n, bw_5_n, bw_6_n, bw_7_n, bw_8_n, bw_12_n, bw_16_n, bw_21_n), names_to = "bw_time_n", values_to = "wt_n")

dat.bw <- dat.clean %>% 
  pivot_longer(., cols = c(bw_0, bw_1, bw_2, bw_3, bw_4, bw_5, bw_6, bw_7, bw_8, bw_12, bw_16, bw_21), names_to = "bw_time", values_to = "wt")

# Create column for day of sampling (numeric)
dat.bwn <- dat.bwn %>% mutate("day" = case_when(bw_time_n == "bw_0_n" ~ 0,
                                                    bw_time_n == "bw_1_n" ~ 1,
                                                    bw_time_n == "bw_2_n" ~ 2,
                                                    bw_time_n == "bw_3_n" ~ 3,
                                                    bw_time_n == "bw_4_n" ~ 4,
                                                    bw_time_n == "bw_5_n" ~ 5,
                                                    bw_time_n == "bw_6_n" ~ 6,
                                                    bw_time_n == "bw_7_n" ~ 7,
                                                    bw_time_n == "bw_8_n" ~ 8,
                                                    bw_time_n == "bw_12_n" ~ 12,
                                                    bw_time_n == "bw_16_n" ~ 16,
                                                    bw_time_n == "bw_21_n" ~ 21))

# Create column for day of sampling (numeric)
dat.bw <- dat.bw %>% mutate("day" = case_when(bw_time == "bw_0" ~ 0,
                                                    bw_time == "bw_1" ~ 1,
                                                    bw_time == "bw_2" ~ 2,
                                                    bw_time == "bw_3" ~ 3,
                                                    bw_time == "bw_4" ~ 4,
                                                    bw_time == "bw_5" ~ 5,
                                                    bw_time == "bw_6" ~ 6,
                                                    bw_time == "bw_7" ~ 7,
                                                    bw_time == "bw_8" ~ 8,
                                                    bw_time == "bw_12" ~ 12,
                                                    bw_time == "bw_16" ~ 16,
                                                    bw_time == "bw_21" ~ 21))
# Order dataframe for analysis
dat.bwn <- dat.bwn[order(dat.bw$day),]

# Remove rows with NA
dat.bwn <- subset(dat.bwn, !is.na(wt_n))
dat.bw <- subset(dat.bw, !is.na(wt))

# Calculate regression lines per treatment group
reg_lines <- dat.bw %>%
  group_by(feedtreat) %>%
  summarize(slope = coef(lm(wt ~ day))[2],
            intercept = coef(lm(wt ~ day))[1])

# Calculate regression lines per rat over time
rl_rats <- dat.bw %>%
  group_by(rat_name) %>%
  summarize(slope = coef(lm(wt ~ day))[2],
            intercept = coef(lm(wt ~ day))[1]) %>%
  cbind(data.frame(feedtreat = rep(c("LF_CTRL","LF_PFOS","HF_CTRL","HF_PFOS"), each = 12),
                   feed = rep(c("LF","HF"), each = 24),
                   treatment = rep(c("CTRL","PFOS"), each = 12, times = 2),
                   day = rep(c("d8","d21"), each = 6, times = 4)))

# Build the linear model
model  <- lm(slope ~ feedtreat, data = rl_rats)
# Create a QQ plot of residuals
ggqqplot(residuals(model))

# Compute Shapiro-Wilk test of normality
shapiro_test(residuals(model))

## # A tibble: 1 × 3
##   variable         statistic p.value
##   <chr>                <dbl>   <dbl>
## 1 residuals(model)     0.952  0.0472

# Test for equal variance
rl_rats %>% levene_test(slope ~ feedtreat)

## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.

## # A tibble: 1 × 4
##     df1   df2 statistic     p
##   <int> <int>     <dbl> <dbl>
## 1     3    44     0.220 0.882

# Save result
EQUAL.VAR <- rl_rats %>% levene_test(slope ~ feedtreat) %>% pull(p) > 0.05

## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.

# Kruskal-Wallis test
pwc1 <- rl_rats %>%
  kruskal_test(slope ~ feedtreat)
pwc1

## # A tibble: 1 × 6
##   .y.       n statistic    df      p method        
## * <chr> <int>     <dbl> <int>  <dbl> <chr>         
## 1 slope    48      7.49     3 0.0579 Kruskal-Wallis

# Posthoc Dunn's test
pwc2 <- rl_rats %>% 
  dunn_test(slope ~ feedtreat, p.adjust.method = "fdr") 
pwc2

## # A tibble: 6 × 9
##   .y.   group1  group2     n1    n2 statistic      p p.adj p.adj.signif
## * <chr> <chr>   <chr>   <int> <int>     <dbl>  <dbl> <dbl> <chr>       
## 1 slope HF_CTRL HF_PFOS    12    12    -0.190 0.850  0.850 ns          
## 2 slope HF_CTRL LF_CTRL    12    12     2.01  0.0442 0.133 ns          
## 3 slope HF_CTRL LF_PFOS    12    12     1.62  0.106  0.158 ns          
## 4 slope HF_PFOS LF_CTRL    12    12     2.20  0.0277 0.133 ns          
## 5 slope HF_PFOS LF_PFOS    12    12     1.81  0.0706 0.141 ns          
## 6 slope LF_CTRL LF_PFOS    12    12    -0.394 0.694  0.833 ns

5.4.2 Statistical analysis

Based on the above tests of the regression for four groups we can analyse the data using Kruskal-Wallis with posthoc Dunn’s test.

# Set variables for inner and outer analysis, and variable for facet
INNER.VAR <- "treatment"
OUTER.VAR <- "feed"
FACETVAR <- "day"

# Statistics for overall data
stat.in2 <- rl_rats %>%
  group_by(.data[[OUTER.VAR]]) %>%
  wilcox_test(as.formula(paste("slope ~",INNER.VAR, sep = " "))) %>%
  adjust_pvalue(method = "BH") %>%
  add_significance("p.adj") %>% 
  p_format("p.adj", accuracy = 0.0001, trailing.zero = TRUE, new.col = TRUE)
stat.in2

## # A tibble: 2 × 11
##   feed  .y.   group1 group2    n1    n2 statistic     p p.adj p.adj.signif
## * <chr> <chr> <chr>  <chr>  <int> <int>     <dbl> <dbl> <dbl> <chr>       
## 1 HF    slope CTRL   PFOS      12    12        73 0.977 0.977 ns          
## 2 LF    slope CTRL   PFOS      12    12        79 0.713 0.977 ns          
## # ℹ 1 more variable: p.adj.format <chr>

stat.out2 <- rl_rats %>%
  wilcox_test(as.formula(paste("slope ~",OUTER.VAR, sep = " "))) %>%
  adjust_pvalue(method = "BH") %>%
  add_significance("p.adj") %>%
  p_format("p.adj", accuracy = 0.0001, trailing.zero = TRUE, new.col = TRUE)
stat.out2

## # A tibble: 1 × 10
##   .y.   group1 group2    n1    n2 statistic       p   p.adj p.adj.signif
## * <chr> <chr>  <chr>  <int> <int>     <dbl>   <dbl>   <dbl> <chr>       
## 1 slope HF     LF        24    24       157 0.00633 0.00633 **          
## # ℹ 1 more variable: p.adj.format <chr>

## Calculate positions statistics on plot
stat.in2 <- stat.in2 %>% add_xy_position(x = OUTER.VAR, dodge = 0.8)
stat.out2 <- stat.out2 %>% add_xy_position(x = OUTER.VAR, dodge = 0.9)
stat.out2$y.position <- max(stat.in2$y.position)*1.1

# Sort dat.clean
rl_rats <- rl_rats[order(rl_rats$feedtreat),]

# Visualise growth rates per group with statistics
p <- rl_rats %>%
  ggboxplot(x = "feed",
            y = "slope",
            color = "feedtreat",
            palette = params$COL1)
p2 <- p + stat_pvalue_manual(stat.in2, label = "p.adj.format",tip.length = 0, hide.ns = FALSE) +
  stat_pvalue_manual(stat.out2, label = "p.adj.format",tip.length = 0, hide.ns = FALSE) +
  scale_y_continuous(expand = expansion(mult = c(0,0.1)))
p2

## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's fill values.

# Output plot
ggsave(filename = paste0("plots/animal_data/weight/bodyweight_slope_combined_plot.png"), p2, device = "png", dpi = 300, units = "mm", width = 100, height = 100)

## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's fill values.

ggsave(filename = paste0("plots/animal_data/weight/bodyweight_slope_combined_plot.pdf"), p2, device = "pdf", dpi = 300, units = "mm", width = 100, height = 100)

## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's fill values.

# Statistics for facet by day
stat.in1 <- rl_rats %>%
  group_by(.data[[OUTER.VAR]],.data[[FACETVAR]]) %>%
  wilcox_test(as.formula(paste("slope ~",INNER.VAR, sep = " "))) %>%
  adjust_pvalue(method = "BH") %>%
  add_significance("p.adj") %>% 
  p_format("p.adj", accuracy = 0.0001, trailing.zero = TRUE, new.col = TRUE)
stat.in1

## # A tibble: 4 × 12
##   feed  day   .y.   group1 group2    n1    n2 statistic      p p.adj
## * <chr> <chr> <chr> <chr>  <chr>  <int> <int>     <dbl>  <dbl> <dbl>
## 1 HF    d21   slope CTRL   PFOS       6     6         9 0.18   0.36 
## 2 HF    d8    slope CTRL   PFOS       6     6        24 0.394  0.394
## 3 LF    d21   slope CTRL   PFOS       6     6        29 0.0931 0.36 
## 4 LF    d8    slope CTRL   PFOS       6     6        12 0.394  0.394
## # ℹ 2 more variables: p.adj.signif <chr>, p.adj.format <chr>

stat.out1 <- rl_rats %>%
  group_by(.data[[FACETVAR]]) %>%
  wilcox_test(as.formula(paste("slope ~",OUTER.VAR, sep = " "))) %>%
  adjust_pvalue(method = "BH") %>%
  add_significance("p.adj") %>%
  p_format("p.adj", accuracy = 0.0001, trailing.zero = TRUE, new.col = TRUE)
stat.out1

## # A tibble: 2 × 11
##   day   .y.   group1 group2    n1    n2 statistic        p    p.adj p.adj.signif
## * <chr> <chr> <chr>  <chr>  <int> <int>     <dbl>    <dbl>    <dbl> <chr>       
## 1 d21   slope HF     LF        12    12        14 0.000371 0.000742 ***         
## 2 d8    slope HF     LF        12    12        55 0.347    0.347    ns          
## # ℹ 1 more variable: p.adj.format <chr>

## Calculate positions statistics on plot
stat.in1 <- stat.in1 %>% add_xy_position(x = OUTER.VAR, dodge = 0.8)
stat.out1 <- stat.out1 %>% add_xy_position(x = OUTER.VAR, dodge = 0.9)
stat.out1$y.position <- max(stat.in1$y.position)*1.1

# Sort dat.clean
rl_rats <- rl_rats[order(rl_rats$feedtreat),]

# Visualise growth rates per group with statistics
p <- rl_rats %>%
  ggboxplot(x = "feed",
            y = "slope",
            color = "feedtreat",
            palette = params$COL1,
            facet.by = "day")
p1 <- p + stat_pvalue_manual(stat.in1, label = "p.adj.format",tip.length = 0, hide.ns = FALSE) +
  stat_pvalue_manual(stat.out1, label = "p.adj.format",tip.length = 0, hide.ns = FALSE) +
  scale_y_continuous(expand = expansion(mult = c(0,0.1)))
p1

## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's fill values.

# Output plot
ggsave(filename = paste0("plots/animal_data/weight/bodyweight_slope_",FACETVAR,"_plot.png"), p2, device = "png", dpi = 300, units = "mm", width = 100, height = 100)

## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's fill values.

ggsave(filename = paste0("plots/animal_data/weight/bodyweight_slope_",FACETVAR,"_plot.pdf"), p2, device = "pdf", dpi = 300, units = "mm", width = 100, height = 100)

## Warning: No shared levels found between `names(values)` of the manual scale and the
## data's fill values.

5.4.3 Create figure

# Mean and SD for growth rates for plotting
rl_rats_msd <- rl_rats %>% 
  group_by(feedtreat) %>% get_summary_stats(!!sym("slope"), type = "mean_sd")

# Slope texts for plot
slope_text <- data.frame(label = c(paste0("Rate = ",round(rl_rats_msd$mean[1],2),"(",round(rl_rats_msd$sd[1],2),") g/day"),
                                   paste0("Rate = ",round(rl_rats_msd$mean[2],2),"(",round(rl_rats_msd$sd[2],2),") g/day"),
                                   paste0("Rate = ",round(rl_rats_msd$mean[3],2),"(",round(rl_rats_msd$sd[3],2),") g/day"),
                                   paste0("Rate = ",round(rl_rats_msd$mean[4],2),"(",round(rl_rats_msd$sd[4],2),") g/day")),
                         feedtreat = names(table(dat.bw$feedtreat)),
                         x = c(rep(15, times = 4)), y = c(rep(1.15, times = 4)))

# Summary samples in groups
dat.bwn %>% group_by(across(all_of(PREDICTOR))) %>% get_summary_stats(!!sym(OUTCOME), type = "mean_sd")

## # A tibble: 48 × 7
##    feed  treatment bw_time_n variable     n  mean    sd
##    <chr> <chr>     <chr>     <fct>    <dbl> <dbl> <dbl>
##  1 HF    CTRL      bw_0_n    wt_n        12  1    0    
##  2 LF    CTRL      bw_0_n    wt_n        12  1    0    
##  3 HF    PFOS      bw_0_n    wt_n        12  1    0    
##  4 LF    PFOS      bw_0_n    wt_n        12  1    0    
##  5 HF    CTRL      bw_12_n   wt_n         6  1.30 0.051
##  6 LF    CTRL      bw_12_n   wt_n         6  1.35 0.043
##  7 HF    PFOS      bw_12_n   wt_n         6  1.31 0.029
##  8 LF    PFOS      bw_12_n   wt_n         6  1.32 0.02 
##  9 HF    CTRL      bw_16_n   wt_n         6  1.38 0.045
## 10 LF    CTRL      bw_16_n   wt_n         6  1.48 0.047
## # ℹ 38 more rows

# Setup rectangles
rec.df <- data.frame(x0 = -0.5, x1 = 0.5, x2 = 7.5, x3 = 21.5, y1 = 1.58, y2 = 1.6809, c0 = "a", c1 = "b", c2 = "c", t0 = "T0", t1 = "Dosing period", t2 = "Wash-out period")

# Generate plot
p <- ggplot(data = dat.bwn, aes(x = day, y = wt_n, fill = feedtreat)) +
  geom_rect(data = rec.df, mapping = aes(xmin = x0, xmax = x1, ymin = y1, ymax = y2, fill = c0), fill = "#f0f0f0", inherit.aes = FALSE) + 
  geom_text(data = rec.df, aes(x = x0+(x1-x0)/2, y = y1+(y2-y1)/2, label = t0), size = 2.5, inherit.aes = FALSE) +
  geom_rect(data = rec.df, mapping = aes(xmin = x1, xmax = x2, ymin = y1, ymax = y2, fill = c1), fill = "#ffffa0", inherit.aes = FALSE) +
  geom_text(data = rec.df, aes(x = x1+(x2-x1)/2, y = y1+(y2-y1)/2, label = t1), size = 2.5, inherit.aes = FALSE) +
  geom_rect(data = rec.df, mapping = aes(xmin = x2, xmax = x3, ymin = y1, ymax = y2, fill = c2), fill = "#ffd0d0", inherit.aes = FALSE) +
  geom_text(data = rec.df, aes(x = x2+(x3-x2)/2, y = y1+(y2-y1)/2, label = t2), size = 2.5, inherit.aes = FALSE) +
  geom_vline(data = NULL, xintercept = 8, linetype = "dashed", color = "#909090") +
  geom_vline(data = NULL, xintercept = 21, linetype = "dashed", color = "#909090") +
  geom_line(aes(group = rat_name, color = feedtreat), stat = "identity", alpha = 0.6) +
  geom_smooth(aes(group = feedtreat, color = feedtreat), color = "#090909", method = "lm", se = FALSE, size = 0.5, fullrange = TRUE, linetype = "dashed") +
  geom_boxplot(aes(group = day), outlier.shape = NA) + 
  geom_point(size = 0.5) +
  scale_fill_manual(values = params$COL1, labels = c("HF-CTRL","HF-PFOS","LF-CTRL","LF-PFOS"), name = "Groups") +
  scale_color_manual(values = params$COL1) + 
  scale_y_continuous(name = "Body weight ratio", limits = c(0.9,1.7), breaks = seq(1,1.7,0.2)) +
  scale_x_continuous(name = "Day", breaks = c(0,1,2,3,4,5,6,7,8,12,16,21)) +
  theme_pubr() +
  guides(color = "none", alpha = "none") +
  facet_wrap(.~feedtreat, labeller = labeller(feedtreat = c("LF_CTRL" = "LF-CTRL","LF_PFOS" = "LF-PFOS","HF_CTRL" = "HF-CTRL","HF_PFOS" = "HF-PFOS"))) +
  geom_text(data = slope_text, aes(x = x, y = y, label = label), position = position_dodge(width = 0.9), size = 3.5)

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

p

## `geom_smooth()` using formula = 'y ~ x'

# Save plot
suppressMessages(ggsave(filename = "plots/animal_data/weight/bodyweight_norm_over_time.png", plot = p, device = "png", dpi = 300, units = "mm", height = 150, width = 260))
suppressMessages(ggsave(filename = "plots/animal_data/weight/bodyweight_norm_over_time.pdf", plot = p, device = "pdf", dpi = 300, units = "mm", height = 150, width = 260))

# clear the environment and release memory
rm(list = ls(all.names = TRUE)) #will clear all objects includes hidden objects.
invisible(gc()) #free up memory and report the memory usage.

5.5.1 Prepare data

This section sets the variables to be used and prepares the data if necessary.

# load data 
load("R_objects/animal_data.Rdata")
params <- readRDS("R_objects/animal_params.RDS")
dat.clean <- subset(dat, !is.na(cecum_norm))

# Set names of variables
PREDICTOR <- c("treatment","feed","dissection")
OUTCOME <- "cecum_wtbw"
SUBJECT <- "rat_name"

# Create formula
PREDICTOR.F <- ifelse(length(PREDICTOR) > 1, paste(PREDICTOR, collapse = "*"), PREDICTOR)
FORMULA <- as.formula(paste(OUTCOME,PREDICTOR.F, sep = " ~ "))

# Summary samples in groups
dat.clean %>% group_by(across(all_of(PREDICTOR))) %>% get_summary_stats(!!sym(OUTCOME), type = "mean_sd")

## # A tibble: 8 × 7
##   dissection treatment feed  variable       n  mean    sd
##   <chr>      <chr>     <chr> <fct>      <dbl> <dbl> <dbl>
## 1 d08        CTRL      HF    cecum_wtbw     6 19.6   4.34
## 2 d21        CTRL      HF    cecum_wtbw     6 17.5   3.31
## 3 d08        CTRL      LF    cecum_wtbw     6  6.82  2.28
## 4 d21        CTRL      LF    cecum_wtbw     6  6.85  1.30
## 5 d08        PFOS      HF    cecum_wtbw     6 19.4   3.73
## 6 d21        PFOS      HF    cecum_wtbw     6 18.9   2.39
## 7 d08        PFOS      LF    cecum_wtbw     6  6.96  1.09
## 8 d21        PFOS      LF    cecum_wtbw     6  7.02  1.29

# SUmmary of original weight in grams
dat.clean %>% group_by(across(all_of(PREDICTOR))) %>% get_summary_stats(cecum_wt, type = "mean_sd")

## # A tibble: 8 × 7
##   dissection treatment feed  variable     n  mean    sd
##   <chr>      <chr>     <chr> <fct>    <dbl> <dbl> <dbl>
## 1 d08        CTRL      HF    cecum_wt     6  6.90 1.61 
## 2 d21        CTRL      HF    cecum_wt     6  6.90 1.4  
## 3 d08        CTRL      LF    cecum_wt     6  2.28 0.709
## 4 d21        CTRL      LF    cecum_wt     6  3.1  0.476
## 5 d08        PFOS      HF    cecum_wt     6  6.62 1.27 
## 6 d21        PFOS      HF    cecum_wt     6  7.84 1.37 
## 7 d08        PFOS      LF    cecum_wt     6  2.47 0.375
## 8 d21        PFOS      LF    cecum_wt     6  2.93 0.457

dat.clean %>% group_by(across(all_of(PREDICTOR))) %>% get_summary_stats(cecum_wt, type = "full")

## # A tibble: 8 × 16
##   dissection treatment feed  variable     n   min   max median    q1    q3   iqr
##   <chr>      <chr>     <chr> <fct>    <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl> <dbl>
## 1 d08        CTRL      HF    cecum_wt     6  5.15  9.51   6.55  5.79  7.69 1.89 
## 2 d21        CTRL      HF    cecum_wt     6  5.57  9.03   6.47  5.83  7.80 1.97 
## 3 d08        CTRL      LF    cecum_wt     6  1.16  3.12   2.45  1.91  2.69 0.775
## 4 d21        CTRL      LF    cecum_wt     6  2.42  3.50   3.32  2.74  3.45 0.707
## 5 d08        PFOS      HF    cecum_wt     6  4.72  8.56   6.60  6.18  7.06 0.884
## 6 d21        PFOS      HF    cecum_wt     6  6.41 10.1    7.4   6.98  8.54 1.56 
## 7 d08        PFOS      LF    cecum_wt     6  1.99  3.11   2.43  2.30  2.56 0.268
## 8 d21        PFOS      LF    cecum_wt     6  2.57  3.59   2.72  2.58  3.28 0.69 
## # ℹ 5 more variables: mad <dbl>, mean <dbl>, sd <dbl>, se <dbl>, ci <dbl>

# Create plot
bxp <- dat.clean %>%
  ggboxplot(x = if_else(length(PREDICTOR) > 1, PREDICTOR[2],PREDICTOR[1]),
            y = OUTCOME,
            color = PREDICTOR[1],
            facet.by = if(length(PREDICTOR) == 3) PREDICTOR[3],
            palette = "jco")
bxp

# Test for outliers
dat.clean %>% 
  group_by(across(all_of(PREDICTOR))) %>% 
  identify_outliers(!!sym(OUTCOME))

##  [1] dissection        treatment         feed              rat_org          
##  [5] rat_name          cage              feedtreat         feedtreatday     
##  [9] bw_minus18        bw_0              bw_1              bw_2             
## [13] bw_3              bw_4              bw_5              bw_6             
## [17] bw_7              bw_8              bw_12             bw_16            
## [21] bw_21             bw_gain08         bw_gain021        bw_gain821       
## [25] bloodvol_8        bloodvol_16       bloodvol_21       cecum_wt         
## [29] cecum_wtbw        cecum_norm        liver_wt          liver_wtbw       
## [33] liver_norm        brain_wt_estimate transit_0         transit_7        
## [37] transit_20        dose_total_ml     pfos_total_mg     pfos_liver_ugg   
## [41] pfos_liver_mg     pfos_liver8_pct   pfos_brain_ugg    pfos_brain_uggbw 
## [45] pfos_brain_ug     pfos_brain_pct    pfos_serum8_ugml  pfos_serum8_mg   
## [49] pfos_serum8_pct   pfos_serum16_ugml pfos_serum16_mg   pfos_serum16_pct 
## [53] pfos_serum21_ugml pfos_serum21_mg   pfos_serum21_pct  pfos_feces8_ugg  
## [57] pfos_feces12_ugg  pfos_feces16_ugg  pfos_feces21_ugg  pfos_cecum_ugg   
## [61] pfos_urine8_ugml  pfos_urine16_ugml pfos_urine21_ugml pH_je_up         
## [65] pH_je_down        pH_ileum          pH_cecum          acetic           
## [69] formic            propanoic         m2_propanoic      butanoic         
## [73] m3_butanoic       pentanoic         m4_pentanoic      hexanoic         
## [77] heptanoic         is.outlier        is.extreme       
## <0 rækker> (eller 0-længde row.names)

# Build the linear model
model  <- lm(FORMULA, data = dat.clean)
# Create a QQ plot of residuals
ggqqplot(residuals(model))

# Compute Shapiro-Wilk test of normality
shapiro_test(residuals(model))

## # A tibble: 1 × 3
##   variable         statistic p.value
##   <chr>                <dbl>   <dbl>
## 1 residuals(model)     0.983   0.719

plot(model, 1)

dat.clean %>% levene_test(FORMULA)

## # A tibble: 1 × 4
##     df1   df2 statistic     p
##   <int> <int>     <dbl> <dbl>
## 1     7    40      1.73 0.130

# Save result
EQUAL.VAR <- dat.clean %>% levene_test(FORMULA) %>% pull(p) > 0.05

This all shows that cecum weight data has three none-critical outliers, is normally distribution and has equal variance. Therefore we use a one-way ANOVA test with Tukey’s honest significance test to test data for all factors of interest, while using unpaired t-test for pairwise comparisons.

5.5.2 ANOVA One-Way test

Perform test

If we had equality of variance we can now run a one-way ANOVA tests anova_test() (if we have equal variance) or a welch_anova_test() (if variance vary).

if(EQUAL.VAR) {
  res.aov <- dat.clean %>% anova_test(FORMULA)
  res.aov
} else {
  res.aov <- dat.clean %>% welch_anova_test(FORMULA)
  res.aov
}

## ANOVA Table (type II tests)
## 
##                      Effect DFn DFd       F        p p<.05   ges
## 1                 treatment   1  40   0.258 6.14e-01       0.006
## 2                      feed   1  40 232.040 3.02e-18     * 0.853
## 3                dissection   1  40   0.620 4.36e-01       0.015
## 4            treatment:feed   1  40   0.100 7.54e-01       0.002
## 5      treatment:dissection   1  40   0.256 6.15e-01       0.006
## 6           feed:dissection   1  40   0.713 4.03e-01       0.018
## 7 treatment:feed:dissection   1  40   0.235 6.30e-01       0.006

ANOVA shows that feed has the only significant effect on the overall data.

Perform posthoc test

A significant one-way ANOVA is generally followed up by Tukey post-hoc tests to perform multiple pairwise comparisons between groups. When running relaxed Welch one-way test, the Games-Howell post hoc test or pairwise t-tests (with no assumption of equal variances) can be used to compare all possible combinations of group differences.

if(EQUAL.VAR) {
  pwc <- dat.clean %>% tukey_hsd(FORMULA)
  pwc
} else {
  pwc <- dat.clean %>% games_howell_test(FORMULA)
  pwc
}

## # A tibble: 49 × 9
##    term            group1 group2 null.value estimate conf.low conf.high    p.adj
##  * <chr>           <chr>  <chr>       <dbl>    <dbl>    <dbl>     <dbl>    <dbl>
##  1 treatment       CTRL   PFOS            0  3.99e-1    -1.19     1.98  6.14e- 1
##  2 feed            HF     LF              0 -1.20e+1   -13.5    -10.4   4.64e-13
##  3 dissection      d08    d21             0 -6.18e-1    -2.20     0.969 4.36e- 1
##  4 treatment:feed  CTRL:… PFOS:…          0  6.47e-1    -2.33     3.62  9.37e- 1
##  5 treatment:feed  CTRL:… CTRL:…          0 -1.17e+1   -14.7     -8.73  2.91e-12
##  6 treatment:feed  CTRL:… PFOS:…          0 -1.16e+1   -14.5     -8.58  4.04e-12
##  7 treatment:feed  PFOS:… CTRL:…          0 -1.24e+1   -15.3     -9.38  9.64e-13
##  8 treatment:feed  PFOS:… PFOS:…          0 -1.22e+1   -15.2     -9.23  1.18e-12
##  9 treatment:feed  CTRL:… PFOS:…          0  1.51e-1    -2.82     3.13  9.99e- 1
## 10 treatment:diss… CTRL:… PFOS:…          0  1.15e-3    -2.97     2.98  1   e+ 0
## # ℹ 39 more rows
## # ℹ 1 more variable: p.adj.signif <chr>

5.5.3 Nested analysis & Figure

Testing two variables are used when there is a nested aspect in the analysis, for example difference in treatment at each timepoint, that would have the timepoint as the outer variable and treatment as the inner variable

# Set variables for inner and outer analysis, and variable for facet
INNER.VAR <- "treatment"
OUTER.VAR <- "feed"
FACETVAR <- "dissection"

# Remove samples with incomplete metadata
dat.clean <- dat.clean[!is.na(dat.clean[,INNER.VAR]) & !is.na(dat.clean[,OUTER.VAR]),]

# Set variables as factors
dat.clean[,OUTER.VAR] <- as.factor(dat.clean[,OUTER.VAR])
dat.clean[,INNER.VAR] <- as.factor(dat.clean[,INNER.VAR])
dat.clean[,FACETVAR] <- as.factor(dat.clean[,FACETVAR])

# Cecum_norm
fit <- aov(as.formula(paste(OUTCOME,"~", OUTER.VAR,"*",INNER.VAR, sep = " ")), data = dat.clean)
anova(fit)

## Analysis of Variance Table
## 
## Response: cecum_wtbw
##                Df  Sum Sq Mean Sq  F value Pr(>F)    
## feed            1 1715.39 1715.39 244.1094 <2e-16 ***
## treatment       1    1.91    1.91   0.2713 0.6051    
## feed:treatment  1    0.74    0.74   0.1050 0.7474    
## Residuals      44  309.19    7.03                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

TukeyHSD(fit)

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = as.formula(paste(OUTCOME, "~", OUTER.VAR, "*", INNER.VAR, sep = " ")), data = dat.clean)
## 
## $feed
##            diff       lwr       upr p adj
## LF-HF -11.95613 -13.49838 -10.41389     0
## 
## $treatment
##                diff       lwr      upr     p adj
## PFOS-CTRL 0.3986117 -1.143631 1.940854 0.6050505
## 
## $`feed:treatment`
##                        diff        lwr       upr     p adj
## LF:CTRL-HF:CTRL -11.7081226 -14.597642 -8.818603 0.0000000
## HF:PFOS-HF:CTRL   0.6466221  -2.242897  3.536142 0.9322867
## LF:PFOS-HF:CTRL -11.5575214 -14.447041 -8.668002 0.0000000
## HF:PFOS-LF:CTRL  12.3547447   9.465225 15.244264 0.0000000
## LF:PFOS-LF:CTRL   0.1506012  -2.738918  3.040121 0.9990233
## LF:PFOS-HF:PFOS -12.2041435 -15.093663 -9.314624 0.0000000

# Statistics costumed for facet plotting
## Pairwise comparison for inner variable
stat.in <- dat.clean %>%
  group_by(.data[[OUTER.VAR]], .data[[FACETVAR]]) %>%
  t_test(as.formula(paste(OUTCOME,"~",INNER.VAR, sep = " ")), 
         paired = FALSE, var.equal = EQUAL.VAR, 
         detailed = TRUE, alternative = "two.sided") %>%
  add_significance() %>%
  p_format("p", accuracy = 0.001, trailing.zero = TRUE, new.col = TRUE)
stat.in

## # A tibble: 4 × 19
##   dissection feed  estimate estimate1 estimate2 .y.    group1 group2    n1    n2
## * <fct>      <fct>    <dbl>     <dbl>     <dbl> <chr>  <chr>  <chr>  <int> <int>
## 1 d08        HF       0.132     19.6      19.4  cecum… CTRL   PFOS       6     6
## 2 d21        HF      -1.42      17.5      18.9  cecum… CTRL   PFOS       6     6
## 3 d08        LF      -0.134      6.82      6.96 cecum… CTRL   PFOS       6     6
## 4 d21        LF      -0.167      6.85      7.02 cecum… CTRL   PFOS       6     6
## # ℹ 9 more variables: statistic <dbl>, p <dbl>, df <dbl>, conf.low <dbl>,
## #   conf.high <dbl>, method <chr>, alternative <chr>, p.signif <chr>,
## #   p.format <chr>

## Pairwise comparison for outer variable
stat.out <- dat.clean %>%
  group_by(.data[[FACETVAR]]) %>%
  t_test(as.formula(paste(OUTCOME,"~", OUTER.VAR, sep = " ")), 
         paired = FALSE, var.equal = EQUAL.VAR, 
         detailed = TRUE, alternative = "two.sided") %>%
  add_significance() %>%
  p_format("p", accuracy = 0.001, trailing.zero = TRUE, new.col = TRUE)
stat.out

## # A tibble: 2 × 18
##   dissection estimate estimate1 estimate2 .y.        group1 group2    n1    n2
## * <fct>         <dbl>     <dbl>     <dbl> <chr>      <chr>  <chr>  <int> <int>
## 1 d08            12.6      19.5      6.89 cecum_wtbw HF     LF        12    12
## 2 d21            11.3      18.2      6.94 cecum_wtbw HF     LF        12    12
## # ℹ 9 more variables: statistic <dbl>, p <dbl>, df <dbl>, conf.low <dbl>,
## #   conf.high <dbl>, method <chr>, alternative <chr>, p.signif <chr>,
## #   p.format <chr>

## Calculate positions statistics on plot
stat.in <- stat.in %>% add_xy_position(x = OUTER.VAR, dodge = 0.8)
stat.out <- stat.out %>% add_xy_position(x = OUTER.VAR)
stat.out$y.position <- max(stat.in$y.position)*1.1

# Create plot
p <- ggboxplot(dat.clean, x = OUTER.VAR, y = "cecum_norm",
               color = "feedtreat",
               fill = "feedtreat",
               add = "jitter",
               add.params = list(size = 1),
               panel.labs = list("dissection" = c("Day 8","Day 21")),
               facet.by = "dissection") +
  theme_pubr(legend = "top") +
  scale_fill_manual(breaks = dat.clean$treatment, values = params$COL1) +
  scale_color_manual(breaks = dat.clean$treatment, values = c("HF_CTRL" = "black","HF_PFOS" = "black","LF_CTRL" = "black","LF_PFOS" = "black")) +
  scale_x_discrete(name = "Feed") +
  guides(color = FALSE) +
  theme(axis.title.x = element_blank())

## Warning: The `<scale>` argument of `guides()` cannot be `FALSE`. Use "none" instead as
## of ggplot2 3.3.4.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

p

# Add p-values for inner and outer analysis to plot
p.stat <- p + 
  stat_pvalue_manual(stat.in, label = "p.signif", tip.length = 0, hide.ns = TRUE, color = "red") + 
  stat_pvalue_manual(stat.out, label = "p.signif", tip.length = 0, hide.ns = TRUE, y.position = c(1.45,1.45)) +
  scale_y_continuous(name ="Percentage from mean HF-CTRL", limits = c(0.15,1.5), breaks = seq(0,1.5,0.25), labels = function(x) paste0(x*100-100, "%"))
p.stat

### Create output plot
suppressMessages(ggsave(filename = paste0("plots/animal_data/weight/cecum_weight_nested_",OUTER.VAR,"_",INNER.VAR,".png"), plot = p.stat, device = "png", units = "mm", width = 130, height = 80, dpi = 300))
suppressMessages(ggsave(filename = paste0("plots/animal_data/weight/cecum_weight_nested_",OUTER.VAR,"_",INNER.VAR,".pdf"), plot = p.stat, device = "pdf", units = "mm", width = 130, height = 80, dpi = 300))

# clear the environment and release memory
rm(list = ls(all.names = TRUE)[ls(all.names = TRUE) != "params"])
invisible(gc())

5.6.1 Prepare data

This section sets the variables to be used and prepares the data if necessary.

# load data 
load("R_objects/animal_data.Rdata")
params <- readRDS("R_objects/animal_params.RDS")
dat.clean <- subset(dat, !is.na(cecum_norm))

# Set names of variables
PREDICTOR <- c("treatment","feed","dissection")
OUTCOME <- "liver_wtbw"
SUBJECT <- "rat_name"

# Create formula
PREDICTOR.F <- ifelse(length(PREDICTOR) > 1, paste(PREDICTOR, collapse = "*"), PREDICTOR)
FORMULA <- as.formula(paste(OUTCOME,PREDICTOR.F, sep = " ~ "))

# Summary samples in groups
dat.clean %>% group_by(across(all_of(PREDICTOR))) %>% get_summary_stats(!!sym(OUTCOME), type = "mean_sd")

## # A tibble: 8 × 7
##   dissection treatment feed  variable       n  mean    sd
##   <chr>      <chr>     <chr> <fct>      <dbl> <dbl> <dbl>
## 1 d08        CTRL      HF    liver_wtbw     6  46.9  2.73
## 2 d21        CTRL      HF    liver_wtbw     6  39.3  3.14
## 3 d08        CTRL      LF    liver_wtbw     6  46.3  3.37
## 4 d21        CTRL      LF    liver_wtbw     6  42.2  4.44
## 5 d08        PFOS      HF    liver_wtbw     6  49.1  3.90
## 6 d21        PFOS      HF    liver_wtbw     6  43.4  3.92
## 7 d08        PFOS      LF    liver_wtbw     6  51.6  3.35
## 8 d21        PFOS      LF    liver_wtbw     6  42.2  4.34

# SUmmary of original weight in grams
dat.clean %>% group_by(across(all_of(PREDICTOR))) %>% get_summary_stats(liver_wt, type = "mean_sd")

## # A tibble: 8 × 7
##   dissection treatment feed  variable     n  mean    sd
##   <chr>      <chr>     <chr> <fct>    <dbl> <dbl> <dbl>
## 1 d08        CTRL      HF    liver_wt     6  16.6  2.42
## 2 d21        CTRL      HF    liver_wt     6  15.5  1.68
## 3 d08        CTRL      LF    liver_wt     6  15.7  2.40
## 4 d21        CTRL      LF    liver_wt     6  19.2  2.37
## 5 d08        PFOS      HF    liver_wt     6  16.8  1.99
## 6 d21        PFOS      HF    liver_wt     6  17.9  1.58
## 7 d08        PFOS      LF    liver_wt     6  18.6  3.38
## 8 d21        PFOS      LF    liver_wt     6  17.7  1.88

dat.clean %>% group_by(across(all_of(PREDICTOR))) %>% get_summary_stats(liver_wt, type = "full")

## # A tibble: 8 × 16
##   dissection treatment feed  variable     n   min   max median    q1    q3   iqr
##   <chr>      <chr>     <chr> <fct>    <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl> <dbl>
## 1 d08        CTRL      HF    liver_wt     6  14.3  19.9   15.8  14.9  18.5  3.61
## 2 d21        CTRL      HF    liver_wt     6  13.4  17.3   16.0  14.0  16.6  2.61
## 3 d08        CTRL      LF    liver_wt     6  13.0  18.5   15.6  13.6  17.7  4.11
## 4 d21        CTRL      LF    liver_wt     6  16.9  23.4   18.4  18.0  20.0  2.01
## 5 d08        PFOS      HF    liver_wt     6  13.5  18.9   16.7  16.2  18.3  2.07
## 6 d21        PFOS      HF    liver_wt     6  15.6  20.1   18.3  17.0  18.4  1.38
## 7 d08        PFOS      LF    liver_wt     6  15.1  24.2   17.3  16.6  20.2  3.60
## 8 d21        PFOS      LF    liver_wt     6  15.4  19.8   17.6  16.4  19.4  3.07
## # ℹ 5 more variables: mad <dbl>, mean <dbl>, sd <dbl>, se <dbl>, ci <dbl>

# Create plot
bxp <- dat.clean %>%
  ggboxplot(x = if_else(length(PREDICTOR) > 1, PREDICTOR[2],PREDICTOR[1]),
            y = OUTCOME,
            color = PREDICTOR[1],
            facet.by = if(length(PREDICTOR) == 3) PREDICTOR[3],
            palette = "jco")
bxp

# Test for outliers
dat.clean %>% 
  group_by(across(all_of(PREDICTOR))) %>% 
  identify_outliers(!!sym(OUTCOME))

## # A tibble: 2 × 79
##   dissection treatment feed  rat_org rat_name  cage feedtreat feedtreatday
##   <chr>      <chr>     <chr> <chr>   <chr>    <int> <chr>     <chr>       
## 1 d08        PFOS      HF    R39     R39         20 HF_PFOS   HF_PFOS_d08 
## 2 d08        PFOS      HF    R43     R41         22 HF_PFOS   HF_PFOS_d08 
## # ℹ 71 more variables: bw_minus18 <dbl>, bw_0 <dbl>, bw_1 <int>, bw_2 <int>,
## #   bw_3 <int>, bw_4 <int>, bw_5 <int>, bw_6 <int>, bw_7 <int>, bw_8 <int>,
## #   bw_12 <dbl>, bw_16 <dbl>, bw_21 <dbl>, bw_gain08 <dbl>, bw_gain021 <dbl>,
## #   bw_gain821 <dbl>, bloodvol_8 <dbl>, bloodvol_16 <dbl>, bloodvol_21 <dbl>,
## #   cecum_wt <dbl>, cecum_wtbw <dbl>, cecum_norm <dbl>, liver_wt <dbl>,
## #   liver_wtbw <dbl>, liver_norm <dbl>, brain_wt_estimate <dbl>,
## #   transit_0 <int>, transit_7 <int>, transit_20 <int>, dose_total_ml <dbl>, …

# Build the linear model
model  <- lm(FORMULA, data = dat.clean)
# Create a QQ plot of residuals
ggqqplot(residuals(model))

# Compute Shapiro-Wilk test of normality
shapiro_test(residuals(model))

## # A tibble: 1 × 3
##   variable         statistic p.value
##   <chr>                <dbl>   <dbl>
## 1 residuals(model)     0.980   0.559

plot(model, 1)

dat.clean %>% levene_test(FORMULA)

## # A tibble: 1 × 4
##     df1   df2 statistic     p
##   <int> <int>     <dbl> <dbl>
## 1     7    40     0.225 0.977

# Save result
EQUAL.VAR <- dat.clean %>% levene_test(FORMULA) %>% pull(p) > 0.05

This all shows that normalised cecum weight data has two outliers, is normally distribution and has equal variance. Therefore we use a one-way ANOVA test with Tukey’s honest significance test to test data for all factors of interest, while using unpaired t-test for pairwise comparisons.

5.6.2 ANOVA One-Way test

Perform test

If we had equality of variance we can now run a one-way ANOVA tests anova_test() (if we have equal variance) or a welch_anova_test() (if variance vary).

if(EQUAL.VAR) {
  res.aov <- dat.clean %>% anova_test(FORMULA)
  res.aov
} else {
  res.aov <- dat.clean %>% welch_anova_test(FORMULA)
  res.aov
}

## ANOVA Table (type II tests)
## 
##                      Effect DFn DFd      F        p p<.05      ges
## 1                 treatment   1  40  7.564 9.00e-03     * 1.59e-01
## 2                      feed   1  40  0.762 3.88e-01       1.90e-02
## 3                dissection   1  40 39.215 2.02e-07     * 4.95e-01
## 4            treatment:feed   1  40  0.059 8.09e-01       1.00e-03
## 5      treatment:dissection   1  40  0.627 4.33e-01       1.50e-02
## 6           feed:dissection   1  40  0.001 9.70e-01       3.47e-05
## 7 treatment:feed:dissection   1  40  2.967 9.30e-02       6.90e-02

ANOVA shows that treatment and day of dissection for the samples have significant effects on the overall data.

Perform posthoc test

A significant one-way ANOVA is generally followed up by Tukey post-hoc tests to perform multiple pairwise comparisons between groups. When running relaxed Welch one-way test, the Games-Howell post hoc test or pairwise t-tests (with no assumption of equal variances) can be used to compare all possible combinations of group differences.

if(EQUAL.VAR) {
  pwc <- dat.clean %>% tukey_hsd(FORMULA)
  pwc
} else {
  pwc <- dat.clean %>% games_howell_test(FORMULA)
  pwc
}

## # A tibble: 49 × 9
##    term             group1 group2 null.value estimate conf.low conf.high   p.adj
##  * <chr>            <chr>  <chr>       <dbl>    <dbl>    <dbl>     <dbl>   <dbl>
##  1 treatment        CTRL   PFOS            0    2.93     0.777      5.08 8.9 e-3
##  2 feed             HF     LF              0    0.930   -1.22       3.08 3.88e-1
##  3 dissection       d08    d21             0   -6.67    -8.83      -4.52 2.02e-7
##  4 treatment:feed   CTRL:… PFOS:…          0    3.19    -0.850      7.23 1.65e-1
##  5 treatment:feed   CTRL:… CTRL:…          0    1.19    -2.85       5.23 8.59e-1
##  6 treatment:feed   CTRL:… PFOS:…          0    3.86    -0.178      7.90 6.55e-2
##  7 treatment:feed   PFOS:… CTRL:…          0   -2.00    -6.04       2.04 5.51e-1
##  8 treatment:feed   PFOS:… PFOS:…          0    0.671   -3.37       4.71 9.7 e-1
##  9 treatment:feed   CTRL:… PFOS:…          0    2.67    -1.37       6.71 3.01e-1
## 10 treatment:disse… CTRL:… PFOS:…          0    3.77    -0.265      7.81 7.45e-2
## # ℹ 39 more rows
## # ℹ 1 more variable: p.adj.signif <chr>

Post hoc analysis however show only significant differences between groupings to be between overall treatment, overall dissection days, as well as cross-day comparisons.

5.6.3 Nested analysis & Figure

Testing two variables are used when there is a nested aspect in the analysis, for example difference in treatment at each timepoint, that would have the timepoint as the outer variable and treatment as the inner variable

# Set variables for inner and outer analysis, and variable for facet
INNER.VAR <- "treatment"
OUTER.VAR <- "feed"
FACETVAR <- "dissection"

# Remove samples with incomplete metadata
dat.clean <- dat.clean[!is.na(dat.clean[,INNER.VAR]) & !is.na(dat.clean[,OUTER.VAR]),]

# Set variables as factors
dat.clean[,OUTER.VAR] <- as.factor(dat.clean[,OUTER.VAR])
dat.clean[,INNER.VAR] <- as.factor(dat.clean[,INNER.VAR])
dat.clean[,FACETVAR] <- as.factor(dat.clean[,FACETVAR])

# Cecum_norm
fit <- aov(as.formula(paste(OUTCOME,"~", OUTER.VAR,"*",INNER.VAR, sep = " ")), data = dat.clean)
anova(fit)

## Analysis of Variance Table
## 
## Response: liver_wtbw
##                Df  Sum Sq Mean Sq F value  Pr(>F)  
## feed            1   10.38  10.384  0.4050 0.52780  
## treatment       1  103.04 103.039  4.0190 0.05117 .
## feed:treatment  1    0.80   0.804  0.0314 0.86027  
## Residuals      44 1128.06  25.638                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

TukeyHSD(fit)

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = as.formula(paste(OUTCOME, "~", OUTER.VAR, "*", INNER.VAR, sep = " ")), data = dat.clean)
## 
## $feed
##            diff       lwr      upr     p adj
## LF-HF 0.9302309 -2.015571 3.876032 0.5278019
## 
## $treatment
##               diff         lwr     upr     p adj
## PFOS-CTRL 2.930289 -0.01551291 5.87609 0.0511665
## 
## $`feed:treatment`
##                      diff       lwr      upr     p adj
## LF:CTRL-HF:CTRL 1.1890424 -4.330161 6.708246 0.9389701
## HF:PFOS-HF:CTRL 3.1891000 -2.330103 8.708303 0.4213654
## LF:PFOS-HF:CTRL 3.8605195 -1.658684 9.379723 0.2566549
## HF:PFOS-LF:CTRL 2.0000577 -3.519146 7.519261 0.7684275
## LF:PFOS-LF:CTRL 2.6714771 -2.847726 8.190680 0.5725927
## LF:PFOS-HF:PFOS 0.6714195 -4.847784 6.190623 0.9879985

# Statistics costumed for facet plotting
## Pairwise comparison for inner variable
stat.in <- dat.clean %>%
  group_by(.data[[OUTER.VAR]], .data[[FACETVAR]]) %>%
  t_test(as.formula(paste(OUTCOME,"~",INNER.VAR, sep = " ")), 
         paired = FALSE, var.equal = EQUAL.VAR, 
         detailed = TRUE, alternative = "two.sided") %>%
  add_significance() %>%
  p_format("p", accuracy = 0.001, trailing.zero = TRUE, new.col = TRUE)
stat.in

## # A tibble: 4 × 19
##   dissection feed  estimate estimate1 estimate2 .y.    group1 group2    n1    n2
## * <fct>      <fct>    <dbl>     <dbl>     <dbl> <chr>  <chr>  <chr>  <int> <int>
## 1 d08        HF    -2.20         46.9      49.1 liver… CTRL   PFOS       6     6
## 2 d21        HF    -4.18         39.3      43.4 liver… CTRL   PFOS       6     6
## 3 d08        LF    -5.35         46.3      51.6 liver… CTRL   PFOS       6     6
## 4 d21        LF     0.00726      42.2      42.2 liver… CTRL   PFOS       6     6
## # ℹ 9 more variables: statistic <dbl>, p <dbl>, df <dbl>, conf.low <dbl>,
## #   conf.high <dbl>, method <chr>, alternative <chr>, p.signif <chr>,
## #   p.format <chr>

## Pairwise comparison for outer variable
stat.out <- dat.clean %>%
  group_by(.data[[FACETVAR]]) %>%
  t_test(as.formula(paste(OUTCOME,"~", OUTER.VAR, sep = " ")), 
         paired = FALSE, var.equal = EQUAL.VAR, 
         detailed = TRUE, alternative = "two.sided") %>%
  add_significance() %>%
  p_format("p", accuracy = 0.001, trailing.zero = TRUE, new.col = TRUE)
stat.out

## # A tibble: 2 × 18
##   dissection estimate estimate1 estimate2 .y.        group1 group2    n1    n2
## * <fct>         <dbl>     <dbl>     <dbl> <chr>      <chr>  <chr>  <int> <int>
## 1 d08          -0.970      48.0      49.0 liver_wtbw HF     LF        12    12
## 2 d21          -0.891      41.4      42.2 liver_wtbw HF     LF        12    12
## # ℹ 9 more variables: statistic <dbl>, p <dbl>, df <dbl>, conf.low <dbl>,
## #   conf.high <dbl>, method <chr>, alternative <chr>, p.signif <chr>,
## #   p.format <chr>

## Calculate positions statistics on plot
stat.in <- stat.in %>% add_xy_position(x = OUTER.VAR, dodge = 0.8)
stat.out <- stat.out %>% add_xy_position(x = OUTER.VAR)
stat.out$y.position <- max(stat.in$y.position)*1.1

# Create plot
p <- ggboxplot(dat.clean, x = OUTER.VAR, y = "liver_norm",
               color = "feedtreat",
               fill = "feedtreat",
               add = "jitter",
               add.params = list(size = 1),
               panel.labs = list("dissection" = c("Day 8","Day 21")),
               facet.by = "dissection") +
  theme_pubr(legend = "top") +
  scale_fill_manual(breaks = dat.clean$treatment, values = params$COL1) +
  scale_color_manual(breaks = dat.clean$treatment, values = c("HF_CTRL" = "black","HF_PFOS" = "black","LF_CTRL" = "black","LF_PFOS" = "black")) +
  scale_x_discrete(name = "Feed") +
  guides(color = FALSE) +
  theme(axis.title.x = element_blank())
p

# Add p-values for inner and outer analysis to plot
p.stat <- p + 
  stat_pvalue_manual(stat.in, label = "p.signif", tip.length = 0, hide.ns = TRUE,  color = "red", y.position = 1.28) +
  stat_pvalue_manual(stat.out, label = "p.signif", tip.length = 0, hide.ns = TRUE) +
  scale_y_continuous(name ="Percentage from mean HF-CTRL", expand = expansion(mult = c(0.01, 0.1)), limits = c(0.85,1.3), breaks = seq(0.8,1.3,0.1), labels = function(x) paste0(x*100-100, "%"))
p.stat

### Create output plot
suppressMessages(ggsave(filename = paste0("plots/animal_data/weight/liver_weight_nested_",OUTER.VAR,"_",INNER.VAR,".png"), plot = p.stat, device = "png", units = "mm", width = 130, height = 80, dpi = 300))
suppressMessages(ggsave(filename = paste0("plots/animal_data/weight/liver_weight_nested_",OUTER.VAR,"_",INNER.VAR,".pdf"), plot = p.stat, device = "pdf", units = "mm", width = 130, height = 80, dpi = 300))

# clear the environment and release memory
rm(list = ls(all.names = TRUE)[ls(all.names = TRUE) != "params"])
invisible(gc())